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Abstract 

A coupled-cluster approach for systems of N bosons in external traps is 
developed. In the coupled-cluster approach the exact many-body wavefunc- 
tion is obtained by applying an exponential operator exp{T} to the ground 
configuration 1 0o) - The natural ground configuration for bosons is, of course, 
when all reside in a single orbital. Because of this simple structure of |</>o), the 
appearance of excitation operators T = J2 n =l -^n f° r bosons is much simpler 
than for fermions. We can treat very large numbers of bosons with coupled- 
cluster expansions. In a substantial part of this work, we address the issue 
of size consistency for bosons and enquire whether truncated coupled-cluster 
expansions are size consistent. We show that, in contrast to the familiar sit- 
uation for fermions for which coupled-cluster expansions are size consistent, 
for bosons the answer to this question depends on the choice of ground con- 
figuration. Utilizing the natural ground configuration, working equations for 
the truncated coupled-cluster with T = T\ + T2, i.e., coupled-cluster singles 
doubles (CCSD) are explicitly derived. Finally, an illustrative numerical ex- 
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ample for a condensate with up to N = 10000 bosons in an harmonic trap is 

provided and analyzed. The results are highly promising. 
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I. INTRODUCTION 



Following the experimental demonstrations of Bose-Einstein condensates in dilute gases 
[1,2], the problem of many bosonic atoms interacting in a trap potential has attracted an 
accelerated interest by the scientific community, see [3,4] and references therein. There are 
many phenomena trapped bosons exhibit that can be described quite well by the standard 
mean-field approach, namely Gross-Pitaevskii (GP) theory [5], see [3,4] and reference therein. 
Side-by-side, the necessity to go beyond mean-field and describe many-body facets of trapped 
bosons has become well-accepted and perused by the community, see the review [6] and 
references therein. 

The many-boson problem is difficult to tackle. Consider, for instance, the standard 
configuration-interaction (CI) approach which employs a basis set expansion. When the 
interaction between the N bosons is substantial and/or many of them are present, the 
number of configurations necessary to properly describe the correlated wavefunction quickly 
increases beyond computational reach and truncations become a must. When truncations 
of the CI expansion are made, there are hints and evidences to slow convergence of the 
CI expansion, see, e.g., [7,8]. Evidently, development of other many-body methods which 
truncate the full configuration space in a different manner are of high relevance and actuality. 

Coupled-cluster theory was first formulated in nuclear physics by Coester [9] and Coester 
and Kiimmel [10], and soon after was introduced to electron-structure theory by Cizek [11] 
and Cizek and Paldus [12]. Coupled-cluster theory has since proven to be a very valuable 
and accurate approach in the many-fermion problem, see [13-15] and references therein. For 
atomic and molecular systems, coupled-cluster theory is currently considered to be one of 
the if not the most powerful many-body tool for calculating electron-correlation energies 
[13-15], also in relativistic systems [16]. In the coupled-cluster approach the exact many- 
body wavefunction is obtained by applying an exponential operator exp{T} to the ground 
configuration \<po). In practice, one truncates of course the operator T. For fermions, it is 
widely known that truncated coupled-cluster expansions are size consistent, which is another 
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advantage the coupled-cluster approach possesses in comparison to truncated CI expansions 
which are not size consistent [17]. 

Our aim in this work is to derive a coupled-cluster theory for bosons with emphasis 
on systems of interacting indistinguishable bosons in traps with up to many particles. We 
investigate aspects like size consistency and what to use as the initial ground configuration 
\4>o). We would like to mention that coupled-cluster approaches for molecular vibrations [18], 
"bosonic nuclei" [19], the spin-boson model [20], and within bosonization of many-electron 
systems [21] have been studied in the literature, but are very different from the present work. 

The structure of the paper is as follows. In section II, we briefly discuss the standard 
configuration-interaction approach. In section III, the coupled-cluster theory for bosons is 
developed, where the issue of size consistency is extensively analyzed. Working equations for 
a truncation of the coupled-cluster to single and double excitations (CCSD) are derived in 
section IV, and an illustrative numerical example is provided in section V. Finally, summary 
and conclusions are drawn in section VI. 

II. THE USUAL APPROACH: CONFIGURATION INTERACTION 

Consider a system of interacting N identical particles, for simplicity either spinless or all 
of the same spin projection. We introduce M one-particle functions (fi(r),i = 1, 2, . . . , M, 
which are called orbitals. The N particles can be distributed over these orbitals and each 
allowed distribution defines a configuration ^i 1 ,i 2 ,...,ijv- ^ ^he particles are fermions, the 
configuration is a determinant 

®h,i 2 ,...,i N = A(p h (p i2 . . . (p iN (1) 
and if they are bosons, it is a permanent 

$*!,«,...,** =SViiVi* (2) 

A and <S denote the antisymmetrizing and symmetrizing operators, respectively. In the 
absence of interaction, each configuration is an eigenfunction of the Hamiltonian H of the 
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system. In the presence of interaction between the particles, the exact eigenfunction \1> in the 
space defined by the M orbitals is given by a superposition of all the allowed configurations 



where the D's are complex numbers. The D's are usually determined variationally by 
diagonalizing the Hamiltonian matrix {(&i lt i 2t ...,i N 1-^1 ^j'i,j2,-,j'jv)}- Clearly, to obtain the 
correct exact eigenfunction, the orbital basis should be complete, i.e. M — > 00, but in 
practical calculations M is kept finite. 

How many distinct configurations participate in the configuration interaction (CI) ex- 
pansion (3)? Two fermions cannot reside in a single orbital and, therefore, the number of 
configurations is simply given by 



In the case of bosons there is no restriction on how many particles can reside in an orbital. 
We find that the number of bosonic configuration reads 



These numbers grow rapidly with the size M of the orbital basis and much more rapidly for 
bosons than for fermions. Consider, for example, 165 particles. For fermions M = 165 is 
needed in order to have a single configuration. Adding just 5 more orbitals, i.e. M = 170, 
increases the number of configurations to over a billion (10 9 ). For bosons, M — 1 is needed 
to have a single configuration and employing M = 170 leads to an astronomically large 
number of configurations. 

For 165 fermions to have only 5 additional (so called virtual) orbitals at their disposal is 
usually insufficient for the calculation of their correlation energy. For an accurate calculation 
more virtual orbitals are required making the CI approach impractical. Fortunately, the 
number of orbitals needed for accurate calculations for bosons is much less than for fermions. 
Because many or even all bosons may reside in a single orbital, the structure of the orbitals 



(3) 



h,i2, ...,ijv 




(4) 




(5) 
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used play a major role in the calculation and the appropriate choice of the orbitals is essential. 
The orbitals are preferentially determined self-consistently as done, for instance, by the use 
of the GP equation [5], see also [3,4]. Nevertheless, to achieve meaningful results M is not 
small and the number of configurations is often beyond reach. To return to our example of 
N = 165, the number of configurations exceeds a billion with just M — 6, i.e., with just 5 
additional (virtual) orbitals. Note that the numbers of bosonic and fermionic configurations 
are identical for the same number of virtual orbitals (M — N orbitals for fermions and M — 1 
for bosons) as can be seen from Eqs. (4) and (5). 

In the following we concentrate on bosons and make use of the destruction and creation 
operators and b\, i = 1,2,...,M, corresponding to the orbitals ipi introduced above. 
These operators fulfill the usual commutator relations 

k, $\=6ij [b i ,b j ]=[b\,b]\=0 (6) 

for bosons. Utilizing these operators, we define the ground configuration 

I0o> = ^(&IHO>, (0o|0o) = l (7) 

which is the ground state of the system in the absence of interaction between the particles. |0) 
denotes the vacuum. All other configurations |3>zi,i 2 ,...,ijv) are obtained directly by applying 
excitation operators to |0). Singly excited configurations read b\b\ \(f>o), doubly excited ones 
are given by b\b^(bi) 2 \(p ), and so on. In analogy to Eq. (3) the exact state can be 
expanded in these orthogonal configurations 

M M 

l*> = I0o> + E d iAbi \<h) + E d ili2 bibi(Pif\ct>o) 

il=2 ii,j2=2 
M 

+ •••+ E d tlt2 ... lN bl 1 bi.--bl(b 1 ) N \ ( j )0 ). (8) 

U,J2,...,ijv=2 

For later use we choose explicitly intermediate normalization of the exact state, i.e., 
(0o I ^) — 1j do not impose normalization on the orthogonal configurations except on \4>o), 
and allow for redundancy in that the same configuration may appear several times in the 
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expansion. We note that the expansion coefficients are independent of the order of the in- 
dices: <l j .j...; x = e/ M ... N . Obviously, there is a one to one correspondence between the 
coefficients in Eq. (8) and those in the expansion in distinct normalized configurations. 

With the aid of the expansion (8) it is relatively straightforward to express the exact 
energy E in terms of the expansion coefficients. Starting from the Schrodinger equation 
H \ = E one immediately arrives at 

E = (0o \H\ *) . (9) 

As usual the system's Hamiltonian consists of an one-particle operator h(f) and a two- 
particle interaction V(f — f ). Expressed in destruction and creation operators H takes on 
the common appearance [22] 

H = J2 Kb\b 3 + i]T v m b\b%\n (io) 

where 




V m = JI tf(f)<P*j(f)V(f- n<Pk(7)<Pi(?Wdr*. (11) 
Inserting Eqs. (8) and (10) into (9) leads to 

M M 

[E - (0o 1^1 0o)] /iV = E [ h n + ( N ~ Waul di + (N-l) £ V llkl d kl . (12) 

1=2 k,l=2 

The energy correction per particle due to the dressing O — > ^ can be expressed by the 
coefficients {di} and {c^}. The orbitals can be conveniently chosen to simplify Eq. (12) 
further by eliminating the {di}; see next chapter for details. 

We should keep in mind that in spite of the compactness of expression (12) this equation 
cannot be used to determine the unknown coefficients {di, dki} since (0o \H\ in Eq. (9) is 
not subject to a variational principle. These coefficients are determined by diagonalizing the 
Hamiltonian matrix which is - as discussed above - of immense dimensionality and in general 
not amenable to practical calculations. One, therefore, resorts to approximations like keeping 
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M very small and/or truncating the CI expansion (8). A particularly appealing approach 
is to truncate the expansion by taking into account only a few classes of configurations. 
In analogy to electron structure calculations we may consider \<f> ) and all singly excited 
configurations (CIS), add to these all doubly excited configurations (CISD), and so on. 

III. COUPLED-CLUSTER THEORY FOR BOSONS 
A. General aspects 

The CI approach discussed in the preceding section is formally straightforward but im- 
practical. Truncating the CI expansion cannot be expected to solve the problem satisfactorily 
in many cases. In particular, when the interaction between the bosons is substantial and/or 
many bosons are present, numerous highly excited configurations may contribute rendering 
systematic truncations impossible. We are, therefore, searching for more efficient approaches 
which are amenable to systematic approximations. 

In the coupled-cluster approach the exact wavefunction is obtained by applying an ex- 
ponential operator to the ground configuration (7): 

l^) = e T |0 o ). (13) 
The operator T is a superposition of excitation operators 

N 

T = £ T n , (14) 

n=\ 

where for bosons we may write 

M 

tn= £ C hi2-inb\ 1 bi2--- b in- ( 15 ) 

iij...jin = 2 

For simplicity we have again introduced redundancies to avoid unpleasant restrictions on 
the summation indices. The yet unknown coefficients c^...^ do not depend on the ordering 
of the subscripts, i.e., c ili2 = c i2il etc. It is convenient to note that [T n ,T m ]=0 and hence 
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exp(T) = exp(T/v) . . .exp(Ti). Because of the simple structure of \<f>o), see Eq. (7), the 
appearance of the coupled-cluster operator T for bosons is much simpler than that for 
fermions, see, e.g., [11,13]. 

Using Eq. (13) and the Schrodinger equation it is easily seen that the exact energy reads 



En = 



e~ T He T 



(16) 



The wavefunction (13) is subject to intermediate normalization (0 O | = 1 as can be 
deduced directly from (0 o |exp(±T) = (0 O |- In this respect the situation is similar to that 
discussed in the preceding chapter, see Eq. (9). On the other hand, Eq. (16) is much more 
powerful because exp(— T)H exp(T) can be evaluated using the useful expansion 

A = e~ T Ae T = A + 1 [A,T] + 1 [[A, T] , T] + . . . , (17) 

which can be applied to any operator A. 

As discussed in the preceding chapter, an expression like (16) is not subject to a vari- 
ational principle and cannot be used to determine the unknown coefficients c ili2 ___ in . To 
proceed we notice that e~ T He T \<f> ) = E |0 O ) and hence projecting on any excited config- 
uration provides an equation for the coefficients. The singly excited configurations lead to 
the (M — l)-equations 

(0o \b\bie- T He T \ <f> ) = % = 2, 3, . . . , M (18) 

and the doubly excited ones to the M(M — l)/2 distinct equations 

% {b\) 2 b i b j e- T He T O ) = i > j = 2, 3, . . . , M (19) 



and so on. The number of independent equations corresponds exactly to the number of 
distinct coefficients, M — 1 coefficients q, M(M — l)/2 coefficients c^-, etc. The equations 
above are nonlinear in the unknown coefficients and, furthermore, are coupled to each other. 
The set (18) contains the q as well as the c^-, while the set (19) also the c^. Since the 
highest possible excitation is iV-fold, the final set of equations does not contain new unknown 
coefficients. 



The total number of distinct coefficients in Eq. (15) is, of course, identical to the total 
number of distinct bosonic configurations given in Eq. (5). We have argued above that this 
number is enormous. Moreover, the equations like (18) used to determine the coefficients 
are nonlinear. So where is the gain with respect to the standard CI method discussed in 
the preceding chapter? The gain is in the favorable properties of the coupled-cluster ansatz 
(13) when truncating the sum of excitation operators in Eq. (14). Let us for demonstration 
include only single and double excitation operators in T, i.e., T = Ti + T 2 . Then, the only 
coefficients available are the q and c^- which can be determined from Eqs. (18) and (19). 
By inspecting that 



one readily notices that this expansion of the wavefunction contains all possible distinct 
configurations of the system. The Eqs. (18) and (19) determine the q and c^- coefficients 
such that the expansion (20) is optimal in providing e~ T He T \ <f> ). In contrast to this CCSD 
approach as we would call it in analogy, the C1SD expansion, on the other hand, knows 
only singly and doubly excited configurations, i.e., is rather related to truncating (20) as 
1^) = I0o) + (^i + -^Tl + T 2 ) \4>$). For additional advantages of the coupled-cluster ansatz 
see the following two sections. 



The influence of Ti is particularly transparent. For this purpose we consider exp(Ti) \4> ) 
and remind that the ground configuration is particularly simple for bosons, see Eq. (7). 
Using the series (17), it is easily seen that 



> = |0 O ) + (7\ + T 2 ) |0 O ) + ^ (Tl + 2T\T 2 + Ti) |0 O ) + . . . 



(20) 



B. Impact of the single excitation operator T\ 



M 



e Tl b\e- Tl =&! + $>,&[ 



(21) 



1=2 



which defines a new creation operator 




(22) 
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where |c| 2 = \°i\ 2 -> which fulfills the boson commutator relation [61,61] = 1. Conse- 
quently, the action of exp(Ti) on the ground permanent O is to define a new permanent 

O )^ e Tl |0 o >oc (61)^10) (23) 

which is, however, not normalized to 1, but rather to (^q | 0o) = [1 + M 2 ]^ 

To proceed, one can consider the quantities appearing in Eq. (22) as the first column 
of an unitary matrix U which defines a new set of M creation operators (b\, 62, ... , = 
(b\,b\, ■ ■ ■ ,b\^XJ. This transformation defines a new set of corresponding orthonormal 
orbitals <f±, <p 2 , ■ ■ ■ , <Pm- In turn, the new set of creation and destruction operators or, 
equivalently, of orbitals, can be formally utilized to eliminate T\ from T. The remain- 
ing operators of T, the T n with n > 2, are now defined with the operators b\ and b\, e.g., 

Clearly, the impact of T\ is to introduce a new orbital ipi optimal for the coupled-cluster 
expansion. In particular, if we put all T n = 0,n > 2, this new orbital can be constructed 
explicitly. As discussed in chapter IV, this orbital then minimizes the energy functional 

(0O|#|0O>. 

C. Size consistency 

Let us consider a super system consisting of R noninteracting replica of our original TV- 
particle system. Clearly, the exact energy of this super system is E (R) = R - E , where E 
is the energy of the iV-particle system. This result will, of course, be reproduced if either the 
full configuration interaction expansion (8) or the coupled-cluster expansion (13-15) is used. 
In general, the full expansion cannot be utilized and one has to resort to approximations. We, 
therefore, have to pose the question whether truncated CI and coupled-cluster expansions 
for bosonic systems lead to energies which scale correctly with the number of replica R, i.e., 
whether these truncated expansions are size consistent. 

Size consistency plays as important role in electronic structure calculations [17]. Imagine, 
for instance, a molecule which is being broken up into fragments or a cluster consisting of 
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weakly interacting atoms. The computational methods used must be size consistent in order 
to describe correctly the break up of the molecule into fragments or the cluster. Indeed, it is 
well known that truncated CI expansions are generally not size consistent whereas truncated 
CC expansions are size consistent for electrons. In the following we would like to address 
the issue of size consistency for bosons. The concept of size consistency is also relevant for 
bosons. Bosonic systems, e.g., in an external double-well trap can be fragmented [23,24], and 
the computational method used must be able to describe fragmentation correctly. Another, 
even more extreme example is the superfluid to Mott-insulator transition of Bose-Einstein 
condensates in a lattice trap [25,26]. In the superfluid phase all bosons communicate with 
each other and in the insulator phase each potential well of the lattice contains a single 
boson which hardly interacts with the other bosons. 

The ground configuration O of the super system is a symmetrized product of the R 
ground configurations of the individual replica. We write 



where b\ k is the creation operator for bosons in the occupied orbital of the k-th replica. The 
Hamiltonian of the super system is, of course, just the sum of the individual Hamiltonians 



We first show that the truncated CI expansion is not size consistent. For this purpose we 
proceed in analogy to the considerations done for fermions (electrons) [17] and assume each of 
the R replica to consist of two orbitals, or equivalently two destruction operators b lk and b 2k , 
of different spatial symmetry. It is sufficient to demonstrate that CID is not size consistent. 
This implies that the expansion of the total wavefunction consists of the superposition of 
the ground configuration |0 O ) given above in Eq. (24) and of the doubly excited configurations 
(blk) (^!fc) 2 l^o) , k — 1, 2, . . . , R. In this space the Hamiltonian matrix 7i representation of 
H is an "arrow" matrix of dimension R + l, the elements of which read 




(24) 



H = H 1 + H 2 + ... + H R . 



(25) 
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n 



Ok 



C 



Hh 



Hkk' = C ( 0o 



{b\ h )\b 2k f H k (b\ k )\b lk y 



) $kk', 



(26) 



where C is the normalization constant of a double excited configuration. Note that all 
matrix elements Hot, k — 1, 2, . . . , R are identical to each other and so are all the Hkk- We 
put for convenience Hok = V and Hkk = (0o \H \ 0o) + A. The diagonalization of 7i can be 
performed analytically by searching for the roots of E — (0 O \H\ 0o) = J2k \Hok\ / [E — Hkk}- 
This immediately leads to 



E o ( J R) = (0o|^|0o) + y- J R 1/2 



V 2 + 



AR 



1/2 



(27) 



which implies that the truncated CI expansion is not size consistent. Using Eq. (24) and 
(25) one sees that the expectation value (<p \H\ </> ) is size consistent and the correction term 
Eq — (0o \H\ 0o) scales as R}l 2 for large R instead of being proportional to R. 

In contrast to the truncated CI expansion, the truncated coupled-cluster expansion is 
size consistent. The operator T is a sum of for the k — 1, 2, . . . , R replica. Each of the 
T( fc ) has the appearance as in Eqs. (14) and (15) for the individual replica. One has just to 
index the destruction and annihilation operators appearing there by a further subscript k 
for the k-th replica. The values of the coefficients c in Eq. (15) are, of course, the same for 
all replica. Clearly, the various commute with each other and, consequently, exp{T} 
can be factorized as ElfcLi exp{T"( fe )} leading to 



l*> = 



N~ 



r (2) 



N~ 



T (R) 



a. 



N' 



|0> 



(28) 



which is size consistent for any truncation of the T^ k \ 

In spite of the favorable structure (28) a major problem arises. If we a priori know that 
our system consists of R noninteracting replica, we may, of course, use the 0o m Eq. (24) 
and obtain a size consistent result. However, the intension is to apply the coupled-cluster 
method not knowing a priori how our system behaves, i.e., whether it is superfluid or breaks 
up into weakly interacting subsystems. Lacking this knowledge, we cannot use the ansatz 
(24) for 0o- Resorting to 
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w ^ (65) ™ l0> (29) 

which does not distinguish between the R replica as is the case in Eq. (24), we may again 
pose the question: is a truncated coupled-cluster ansatz size consistent? 

To proceed, we first have to identify the b\ operator appearing in Eq. (29) in terms of 
the operators b\ k of the individual replica. Since all replica are identical, we can construct R 
new operators B\ x , Bl 2 , . . . , B[ R of the super system by linearly combining the b\ k . Without 
loss of generality we can always chose 

b\ ee Bl = RrW (bl +b\ 2 + ... + bQ , (30) 

i.e., as a trivial superposition of the creation operators corresponding to the occupied 
orbitals of the individual replica. All the B\ { will posses different permutational sym- 
metries which simplifies the evaluation considerably. For instance, for R = 2 we have 
B[ x = 2-V 2 (V u + b{ 2 ) and B\ 2 = 2" 1 / 2 (V u - b{ 2 ). We note that for each set of virtual 
orbitals an analogous procedure can be applied to introduce the remaining orbitals of the 
super system: b 2l , b 22 , • • • , b\ R are linearly combined to give B 21 , B\ 2 , . . . , B\ R and so on. 
This results in R • M creation operators of the super system emerging from the M operators 
of each of the replica. Since only one orbital is occupied in the ground configuration of the 
super systems, all the other ones are virtual orbitals, i.e., also the Bi2,Bl 3 , . . . ,B\ R refer 
now to virtual orbitals. 

The Hamiltonian (25) and the coupled-cluster operator T are now expressed in the B\ k 
of the super system. Let us consider as an example the one-body part of H in the occupied 
space of the individual replica: 

R R 

J2 h nb\ k b lk =J2hnB\ k B lk . 

k k 

Note that in the two-body part of the H operator products like B\ k B\ k ,B lk B lk i,k ^ k', 
appear. Let us begin the analysis by inspecting the mean-field energy {(p \H\(f) ). Here, 
only the terms of the Hamiltonian containing B^Bn and B^B^BuBn contribute. These 
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terms take on the explicit appearance 

where hu and Vnn are the quantities defined in Eq. (11) for an individual replica. One 
immediately finds 

(0o \H\ O > = NRh n + NR(N * ~ 1 V im 

^{^^IM^} (31) 

implying that even the mean- field energy is not size consistent; a surprising result. The 
mean-field energy of an individual replica is Nhu + JV ( J ^~ 1 ) y im . Consequently, size consis- 
tency is achieved only if each individual replica contains many bosons, i.e., for N 1. 

To better understand the implications of the above finding, let us briefly consider the 
coupled-cluster operator T of the super system. Since now there is only a single occupied 
orbital (related to -Bn), all the other operators B ik relate to virtual orbitals of the super 
system. Consequently, T can be broken up into a part T' which contains excitations solely 
within the original occupied orbitals of the different replica and the remaining part T" where 
excitations to the originally virtual orbitals of these replica are included. As an example we 
consider the double excitation operator T 2 (see Eqs. (14) and (15)): 

r „ R M 

T 2 = T> + T' 2 ' = (B\ k ) (B n ) 2 + £ Y:^B\ k B] k (B u ) 2 . (32) 

k=2 k=l i,j 

In T' 2 ' the terms with Bn and those of V are not included as indicated by the primed 
summation symbol J2'- m the example of two replica, we have the T 2 excitation operator 
(&ti — &i 2 ) (^ii + which is actually an excitation within the occupied manifold of the 
replica. Obviously T' and T" commute. 

Interestingly, the full impact of exp{T'} is needed in order to restore the size consistency. 
Indeed, a calculation shows that 



e- T 'He r 



R { Nhn + N ( N ^V im \ (33) 
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which is the expected correct mean-field result and is identical to the expectation value of 
H obtained with the ansatz (24) for O where the knowledge of having R replica has been 
used. 

The result (33) follows only if the expansion of exp{T'} is fully considered and not 
truncated. The impact of exp{T'} is to transform O in Eq. (29) into the form of Eq. (24) 
which is appropriate for R replica. In other words, truncated coupled-cluster expansions are 
not size consistent once the ansatz (29) is used for <f> . The good news is that the violation 
of the size consistency at least for the mean-field energy leads to negligible errors for large 
individual systems (N ^> 1), see Eq. (31). In this respect bosons and fermions behave 
differently. Due to the fact that each fermion resides in its own orbital, size consistency in 
truncated coupled-cluster expansions follows straightforwardly, as was also found above for 
bosons starting with the O of Eq. (24). 



D. On the choice of the ground configuration 4>q 

In contrast to fermions, the choice of the structure of the ground configuration O as the 
starting point is crucial for bosons if fragmentation or, in particular, phase transitions like 
the superfluid to Mott-insulator transition are to be studied. In the absence of interaction 
between the bosons, the exact ground state has the appearance oc (b\^j |0). It is, therefore, 
natural to start in the presence of interaction from an analogously structured ground con- 
figuration 0o as done in what follows Eq. (7). In the presence of interparticle interaction we 
have the freedom to choose the orbitals defining the destruction and annihilation operators. 
At least as long as this interaction is weak, it is favorable to choose the occupied orbital 
which minimize the energy expectation value (0o \H\ <f>o). This readily leads to the equation 

h + (N — 1) J n ] Mr) = ^i(f) (34) 

which determines the occupied orbital (pi(f). The number /ii can be called orbital energy 
or chemical potential. The direct interaction operator J n is a local operator and reads 
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J n = J <p\(f)V(f- f>i(0^. (35) 

Eq. (34) defines an hermitian Fock-like operator 

F = h + (N-l)J n 
Ftpi = (36) 

the eigenfunctions of which define a complete set of orthogonal orbitals to be used in the 
coupled-cluster calculation. 

For convenience (see chapter IV) one may introduce the more physical operator F 

Jn + Kn 

Fifi = law (37) 
which also contains the nonlocal exchange interaction operator Kn: 

Knfi = J <ft(f)V(f- f>*(f>i(0^- (38) 

Because of the structure of O! both F and F produce the same occupied orbital tpi and 
the same chemical potential. All other orbitals and orbital energies are generally different. 
To avoid confusion, we shall indicate in the following which set of orbitals has been used. 
Finally, we would like to point out that if one chooses V(f—r / ) oc 5(f—r f ), both F(fi = 
and F<pi = /ii<fi reduce to the well-known and widely used GP equation [3,4]. 

As long as the system does not undergo a break up like in the superfluid to Mott-insulator 
transition in an optical lattice potential, 0o of Eq. (7) and the orbital set of Eq. (36), or, 
preferentially, of Eq. (37) can be used in the coupled-cluster calculations. What to do when 
a break up is possible? Here, we would like to stress that Eq. (34) has been obtained from 
the minimization of the mean-field energy (<f> \H\ <f> ) within the ansatz (7) for O - But, 
this ansatz does not necessarily lead to the lowest possible mean-field energy, i.e., it is not 
necessarily the best mean-field ansatz. The best mean-field ansatz allows the bosons to 
reside in different orbitals [27]: 
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h + 



|0o) oc (b\) nr . . . (4)" 2 (b\) ni |0> , ni +n 2 + ...n r = N. (39) 

The number r of different orbitals as well as the occupation numbers n^i = 1,2, ...,r 
which tell us how many bosons reside in which orbital, are not a priori fixed numbers but 
are determined variationally to minimize the mean-field energy. The r optimal orbitals 
involved are, of course, also determined variationally. For brevity of presentation, we do not 
present the equations of the best mean-field approach and refer to the literature [27]. 

The best mean-field ansatz has been shown to be flexible enough to predict and describe 
fragmentation and superfluid and a whole zoo of insulator phases [23,24,28]. We, therefore, 
have reason to expect that (39) provides a useful starting point for many coupled-cluster 
studies. Other ways to determine the orbitals and their occupation numbers can also be 
anticipated in connection with the coupled-cluster approach. 

IV. DERIVATION OF THE WORKING EQUATIONS 

In this chapter the working equations of the coupled-cluster approach are derived and 
discussed. We concentrate here on the ansatz -±^(b\) N \0) for the ground configuration 
for which all the necessary ingredients have been introduced and discussed in the preceding 
chapter. Working equations can also be derived starting from ansatz (39) for O - This ground 
configuration contains several occupied orbitals and consequently the working equations are 
more elaborate. 

We begin by transforming the boson destruction and creation operators with exp{T}. 
Using the expansion (17) one readily finds that (A = e~ T Ae T ) 

bi = h 

b\ = b\, i = 2,3,...,M. (40) 

The destruction operator corresponding to the orbital occupied in O , and the creation 
operators of the virtual orbitals are invariant to the coupled-cluster transformation. In 
contrast, the respective dual operators change 
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b\ = b\ - & 

h = h + d (41) 



where 

N 



\n-l 



n=l 
N 

A = E^(M n - (42) 

n=l 

The operators t n can be found in Eq. (15) and the operators t$ operate in the virtual space 
and read 

M 

t (i) = V . . . U\ U\ u\ 

*2,«3v:*n=2 

4° = <k. (43) 
In the calculations below it is gratifying to note that the C operators commute 

[d,Cj] = [A, A] = (44) 
and that their action on (0o| from the right is simple: 

(0 o |A = o, (0 O | (&i) m £i = o, 

(0o I d = Ci (0o I &i- (45) 

To proceed, we break up the Hamiltonian (10) into several terms according to the number 
of operators related to the occupied orbital ipi. The transformed one-body part H of the 
Hamiltonian then consists of four terms 

M M M 

Ho = hublh + £ h lk b\b k + E /ifcifttfti + E /iw&t&j (46) 

fc=2 A:=2 fc,i=2 

out of which the second is the most involved one. The transformed two-body operator 
V contains many contributions which can be casted into nine terms which, for ease of 
presentation, are listed in the Appendix. 
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We now calculate the energy E = (^> H o y, see Eq. (16). The first term of H in 
Eq. (46) and that of V in the Appendix contribute because of Eq. (45) only to the mean-field 
energy giving 



(<t> \H\<p ) = N 



hu H — Vim 



(47) 



The only terms contributing to the energy correction Eq — (<f>o \H\ 0o) are the second term 
of H in Eq. (46) and the second and forth terms of V in the Appendix. The final result for 
the exact energy reads 

r m n-1 M ) 

E = (0 O \H\ O ) + iV E [^i* + ( N ~ Wn*] °k + — E Vim (2c w + c fc c,) . (48) 

[fc=2 ^ fc,/=2 J 

Inspection of this expression makes clear from which of the above mentioned contributing 
terms the various matrix elements originate. The appearance of the result (48) is similar to 
that derived by the configuration interaction approach, see Eq. (12). The major difference is 
in the c^q term which is missing in the CI expression (12) and arises due to the contribution 
of the single excitation operator 7\ to the wavefunction, see Eq. (20). 

Until now the orbitals used are arbitrary and have not been specified. If we utilize the 
optimized orbitals arising from the Fock-like operators F and F discussed in section III.D, 
we obtain in both cases the same results for the exact energy: 

N(N — 1) 

E = (0 O \H\ O ) + >- E Vuki (2c fc , + c fcCj ) . (49) 

1 k,l=2 



The other term in Eq. (48) has disappeared due to the fact that yp^ F ipij = 0, see 
Eqs. (36) and (37). In analogy to the notion of electron correlation energy [17] we might 
call the correction E — (0 O \H\ 0o), which is caused by the interparticle interaction beyond 
the mean field, boson correlation energy. 

To determine the coefficients {c^} and {c^} we have to evaluate the series of coupled 
equations (18), (19) and so on as discussed in section III. A. The series consists of iV sets of 
such equations, one set for each type of excitation operator T n , n = 1, 2, . . . , N. In practical 
calculations the expansion T = J] T n is truncated. For instance, if Ti and T 2 are considered 
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and the T n ,n > 3, are put to zero, then only the sets of equations (18) and (19) must be 
considered in order to determine the derived coefficients. In the following we calculate this 
CCSD approach as we may call it. 

Whereas the expression (49) for the energy is invariant to the choice of either F or 
F to define the orbital basis used, the equations determining the coefficients do depend 
on this choice. We have computed these equations for an arbitrary set of orthonormal 
orbitals, but present here only the results obtained with the orbitals of F. Let us begin with 
{^0o b\biH 0o ) = 0. From the nine terms of V shown in the Appendix it is easy to see that 
the fifth, eight and ninth terms do not contribute as they all exhibit two creation operators 
for virtual orbitals. All other terms contribute. Using Eqs. (41-45) we obtain a set of M — 1 
coupled equations (i — 2, 3, . . . , M) 

M ( 1 M 

(Hi - Hi) Ci = J2 h ik {2c ki + c k d) + {N-l)l-J2 ( V Ulk + V Uk i) c k + 

k=2 { 1 k=2 

M M \ 

( V Ukl ~ VuklCi) {2c H + CfcCj) + (N - 2) V ^kl (3c fcH + CfcQi + CjCjk) \. (50) 

k,l=2 k,l=2 ) 

This set of equations is the result of the exact evaluation of (<f>o b\biH O ) — an d thus 
contains the coefficients ctu of T 3 . These coefficients have to be put equal to zero if CCSD 
is to be evaluated. Consulting section III.B we see that Ck ^ implies the introduction of a 
new optimized orbital and we may assume that in CCS and M — > oo this orbital is just the 
eigenfunction ipi of F (or, equivalently, F). 

To complete the CCSD we have to solve also for the set of M(M — l)/2 distinct coupled 



equations resulting from < 0, 



b\) bAH 



4>qJ = 0, see text around Eq. (19). Here, all terms 
of H contribute except of the third term of Hq in Eq. (46). Using the relations (41-45) and 
the expressions of H given in Eq. (46) and in the Appendix, the derivation of the coupled 
equations is lengthy but straightforward. In principle, one could derive diagrammatic rules to 
simplify the procedure in analogy to the situation for fermions [11,12], but this is unnecessary 
for bosonic systems, at least as long as O in Eq. (7) is used. The resulting set of coupled 
equations reads = 2,3, ... , M): 
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M 



2 (2/ii -pLi- fxj) Cij - Viju = ( v ijki + Vjiki) c k - (VnuCj + VjmCi) + 



k=2 



M 



M 



M 



J2 i( v uik + IW) a kj + (i <-> j)] + XI Viiifc/Jfctj + H ^nw 7fci*j - (51) 

k — 2 k ,/ — 2 

M M 

X (^li^Cj + VijfejCj) CfcCj + X ^ijfci (2c fci + CfcCj) . 
fc,Z=2 fc,Z=2 

In contrast to Eq. (50) which contains c k u coefficients arising from T 3 , we have concentrated 
in Eq. (51) on CCSD and put all coupled-cluster operators T n , n > 3, to zero. The quantities 
a, (3 and 7 appearing in Eq. (51) are given by 



akj = c k j(N - 3) - c k Cj 

fikij = 2 [(ckjCi + CkiCj) (3N - 5) + 2cijC k (N - 2) + CiCjC k } 



(52) 



Tfciii = 4 Cfci qy(iV - 2)(N - 3) - 2(A - 2) 



2Qj (2c M + CfcCj) + CfciQCj + CuC k Cj + 



- (2c kl + c fc Ci) (2cjj - CjC,-) . 



It is worth noting that the equations (50) arising from (0 O b\biH 

blYbibjH 



/ = are all homoge- 
) = are inhomoge- 



neous, whereas the equations (51) originating from \ < ' 
neous. The inhomogeneity V^n is due to the fifth term of V given in the Appendix, i.e., 
from the only term which is invariant to the exp(T) transformation. 

Before closing this chapter let us briefly discuss CCS. Here, we have to put in Eq. (50) 
all the Chu = as well as all the c k \ = and disregard the set of equations (51). The 
resulting equations are homogeneous in the c k coefficients and c k = is a proper solution. 
This implies that CCS leads to that mean-field energy which is the minimum of (0o \H\ (po), 
see section III.D. Would we have not used the orbitals of F but rather some set of arbitrary 
orthonormal orbitals, then Eq. (50) will become an inhomogeneous equation and the c k 7^ 0. 



V. ILLUSTRATIVE EXAMPLE 



As an example we apply the CCSD approach to N interacting bosons in an ex- 
ternal trap and restrict the orbital space to two orbitals ipi and ip 2 of different spa- 
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tial symmetry. Consequently, the CCSD (or equivalently the CCD) wavefunction reads 



|^) = exp 



C22 



b ). In other words, the wavefunction depends only on a single 



unknown parameter €22- It is easily seen that 



N/2 



l*> = E 



m=0 



u 22 

m! 



AH (2m)! 



1/2 



|JV - 2m, 2m) , 



(53) 



(N — 2m)! 

where for simplicity we assume A" to be an even number and \m,m') is the normalized 
configuration with m bosons in ipi and m! bosons in <p 2 - 

We remind that in coupled-cluster theory intermediate normalization of the wavefunction 
is used, (0o I ^) — 1, an d define the norm of the wavefunction 



(54) 



Using Eq. (53) it is readily shown that this norm obeys a local "decay" law as a function of 
the parameter C22: 

dM (n 2 ) 



■ M 



(55) 



<^C 2 2 C 2 2 

where (722) is the expectation value of the occupation number of bosons in orbital <f2- Because 
of the different spatial symmetry of <pi and ip 2 , these orbitals are the eigenfunctions of the 
reduced one-particle density matrix (natural orbitals) and the respective eigenvalues are (rii) 
and (n 2 ) with (rii) + (n 2 ) = N. Clearly, 



(n 2 ) = (* 



(56) 



which can be evaluated by using Eq. (53). Analogously, the variance of the occupation 
number of bosons in orbital y?2 can be obtained from the second derivative of the norm 

d 2 Af (n 2 2 ) - (n 2 ) 



dc 2 



22 



r 2 
c 22 



(57) 



where {n 2 ,) = ( ^ 



(b\b 2 



/ I ^). In the absence of interaction between the bosons, 
C22 = and (712) = 0, M = 1. As the interaction strength grows, the value of the coupled- 
cluster coefficient \c 22 \ grows as well and with it the mean number of bosons in the orbital ip2- 
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The quantity \(n 2 ) / C22 1 increases and determines the rate of change of the norm according 
to Eq. (55). 

To be specific, we consider now the widely-used, one-dimensional harmonic trap poten- 
tial, — + ^x 2 , and use the contact interaction V(x — x') = X 5(x — x'), see [3,4] and 
references therein. We would like to examine here the performance of the CCSD approach. 
It should be reemphasized that the CCSD wavefunction contains only a single parameter 
c 22 . To solve for this parameter Eq. (51) can be used which reduces to the simple quadratic 
equation 

cj 2 a (N 2 - 7N + 9) + c 22 [(/i 2 - //O + a(N - 3) + (3} + a/4 = (58) 

where 



a = A J \(fi(x)\ 2 \(p 2 (x)\ 2 dx 

P = y |v^i(^) | 4 ^ + J \ip 2 {x)\ A dx 



Note that ji\ here is the usual chemical potential of the GP equation. The ground state 
energy of the CCSD approach reads 

Eo(CCSD) = Eqp + N(N - l)ac 22 . (59) 

Here, E GP = (0 O \H\ O ) is the usual ground state GP energy 

N- 1 



Eqp = N 



h n + 



-A J \ifi{x)\ A dx 



(60) 



For completeness we would like to compare our CCSD results with those of the CISD. 
The latter wavefunction also contains only one parameter d 22 (see chapter II) and the ex- 
pression for the energy i? (CISD) is identical to that in Eq. (59) if we replace c 22 by <i 22 . The 
CISD wavefunction is, however, a superposition of only the two configurations \N,0) and 
I N — 2, 2). The value of the parameter can be simply obtained by diagonalizing the Hamil- 
tonian H in the space of these two configurations. This leads to the following quadratic 
equation 
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d 2 22 aN(N - 1) - 2d 22 [(/i 2 - /n) + a(N - 3) + 0\ - a/2 = (61) 

for the configuration interaction parameter. 

The results of our numerical example are summarized in Figs. 1-3. In Fig. 1 we test the 
performance of CCSD method in terms of the correlation energy. The correlation energy is 
defined as the difference between the Gross-Pitaevskii energy, E GP , and the exact energy, 
.E^o (exact). The latter is obtained in our model by diagonalizing the many-body Hamiltonian 
within the full configuration-interaction space spanned by \m, m') , m — 0, N, m' — N—m. 
We first calculate the CCSD energy £ (CCSD) using Eqs. (34-36,58-60). How much of 
the exact correlation energy Eqp — E (ex&ct) is captured by CCSD is given in percent by 
^correlation = 100 • ■ We have calculated %£ correlation for N = 100, 1000 and 

10000 for several values of the interaction strength A . The results are plotted in Fig. 1 
versus the couping constant A = X (N — 1), which is the only interaction parameter entering 
the GP energy, see Eq. (60). For comparison, the corresponding values obtained by the CISD 
method were calculated as well. We remind that both methods contain one parameter only, 
c 22 and d 22 respectively. It is seen that the CCSD is remarkably successful in obtaining 
the correlation energy, with absolute error of less than 4% up to a huge coupling constant, 
A = 2 • 10 4 . The quality of the CISD, on the other hand, starts to deteriorate already from 
A = 1 on and saturates at about 50%, see Fig. I. Another result observed in Fig. 1 is that 
the performance of CCSD in terms of %£ , C orreiation improves with increasing N, whereas that 
of CISD worsens. 

Next, let us examine the many-body wavefunction obtained by the CCSD method and 
compare it to the exact one. For this, we first normalize the CCSD wavefunction (53) and 
express it as Ylm=oC 2m \N — 2m, 2m). In Fig. 2 the absolute value of the C 2m coefficients 
(they alternate in sign because c 2 2 is negative) for N = 10000 bosons and A = 100 are plotted. 
Although the coupling constant is large, it is remarkable that the CCSD C 2m coefficients 
almost perfectly match the exact coefficients and it is difficult to distinguish between the 
red and black curves of Fig. 2. Another property of the many-body wavefunction when A is 
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growing is that the tail of the coefficients Ci m is extending further, showing that more and 
more excited configurations contribute to the many-body wavefunction. For comparison, the 
two coefficients of the CISD are also shown, which deviate much from the exact solution, 
see Fig. 2. 

Finally, we examine the capability of the CCSD method to reproduce the exact ground- 
state depletion, i.e., the average number (712) of bosons occupying the orbital if2- As men- 
tioned above, (712) and (rii) = N — (712) are the eigenvalues of the reduced one-body density 
and hence are a very sensitive tool for the quality of the CCSD many-body wavefunction. 

We have calculated (712) for N = 100, 1000 and 10000 for several values of the interaction 
strength A up to the huge value of A = 2 • 10 4 . The results are plotted in Fig. 3 together 
with the exact ones and the corresponding values obtained with CISD. It is seen that the 
CCSD is extremely successful in obtaining the depletion (n 2 ) up to a large coupling constant, 
A = 10 2 . From about this value on, the quality of the CCSD wavefunction in terms of (n 2 ) 
depends on the specific number N of bosons. For N = 100 it is very good for all values 
of A. For N = 10000 at the extreme value A = 2 • 10 4 it predicts 2.5 as much depletion 
as the exact many-body wavefunction gives, namely almost 20 bosons instead of 8 out of 
10000 bosons. We remind that all these results are obtained with a single parameter c 2 2- 
The deviations of (n 2 ) for large iV and larger interaction strength Ao is related to the tail of 
the C 2m distribution. As N and A increase, there are more and more non-negligible CCSD 
coefficients which start to deviate from the exact ones. While this does not lead to an error 
of more than 3% percent in the correlation energy, see Fig. 1, it does influence the more 
sensitive measure of exactness of the wavefunction, (n 2 ). For comparison, we also computed 
the corresponding (n 2 ) values with CISD and plotted the results in Fig. 3. We obtained that 
the values of (n2) for all iV saturates at about 0.18 with increasing A, which is more than 
an order of magnitude smaller than the exact and CCSD results. This near independence of 
(JI2) in CISD from the number of bosons N is a manifestation of the minimal correlations 
embedded in the CISD wavefunction, in contrast to the CCSD wavefunction. 

Summarizing the results depicted in Figs. 1-3, we see that the CCSD for bosons performs 
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remarkably well even for large interaction strengths. Utilization of the ground configuration 
in Eq. (7) is an appropriate choice for the coupled-cluster expansion at least for this example 
(see also the discussion below). 

VI. SUMMARY AND CONCLUSIONS 

The standard configuration-interaction approach rapidly becomes impractical in the 
many-body problem. When the interaction between N bosons is substantial and/or many of 
them are present, the number of configurations necessary to correctly describe the correlated 
wavefunction quickly increases beyond computational reach. In searching for more efficient 
approaches which are amenable to systematic approximations (truncations), we have devel- 
oped in this paper a coupled-cluster theory for systems of bosons in external traps. 

In the coupled-cluster approach the exact wavefunction is obtained by applying an ex- 
ponential operator exp{T} to the ground configuration |0 O ). The ground configuration |0 O ) 
depends, of course, on the particle statistics. While for fermions it is a determinant with 
M = N different orbitals, the situation for bosons is more intricate. Since there is no 
limitation on the number of bosons occupying a certain orbital, there are ample legitimate 
choices for the ground permanent of N interacting bosons over M available orbitals. The 
most natural choice for non-interacting or weakly-interacting bosons is, of course, to let all 
bosons reside in the orbital lowest in energy tpi, |0 O ) = '^ffi("i) N "which is our main 
choice for the coupled-cluster theory presented here. 

Because of the simple structure of |0 O ), the appearance of excitation operators T = 
J2n=i Tn for bosons is much simpler than for fermions. When the simplest truncation T — T\ 
is chosen, namely CCS, the effect of exp{Ti} on \<j> ) is to transform </?i to another orbital 
(p\. exp{Ti} optimizes this orbital by mixing the M available orbitals. This reminds the 
situation encountered for fermions, where the operation of the fermionic T\ transforms the 
ground determinant into another determinant (Thouless theorem [29]). 

In a substantial part of this work we addressed the issue of size consistency for bosons 
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and enquired whether truncated coupled-cluster expansions are size consistent. It turns out 
that the answer to this question depends on the choice of ground configuration (perma- 
nent). Considering R non- interacting replica of the iV-boson system, it has been found that 
truncated coupled-cluster expansions are not size consistent with the simplest choice for the 
i?-replica ground permanent, \<f> ) = ^/^ R y (&i) |0)> already for the mean-field energy 
(4>o \H\ (f>o). This is a surprising result when compared to the case of fermions. Fortunately, 
this violation of size consistency, at least for the mean-field energy, leads to negligible errors 
for large individual systems (N ^> 1). Can size consistency in bosonic systems be fully 
restored, perhaps with another choice of the i?-replica ground permanent? Yes, it has been 
straightforwardly shown that truncated coupled-cluster expansions are size consistent with 
the ground permanent |0 O ) = jj^r/j Uk=i (&i fc ) |0), which "distinguishes" between the R 
replica; also see discussion below. 

Next, we moved to derive working equations of the coupled-cluster approach for the 
natural ground configuration |0 O ) = ^?(^1) JV First, it has been shown that the exact 
correlation energy depends on two kinds of coefficients only: {c k } and {c^} of the single 
and double excitation operators 7\ and T 2 . For a given truncation of the coupled-cluster 
exponential operator exp{T}, it is possible in principle to calculate the correlation energy. 
Here, for the specific truncation of T = Ti + T 2 , i.e., CCSD, general working equations for 
{cfc}, {cki} have been explicitly derived. 

Finally, we tested the performance of the CCSD in a model where an exact solution 
can be computed. We employed an harmonic trap and restricted the orbital space to two 
orbitals of different spatial symmetry. The exact solution is obtained, of course, by diagonal- 
izing the many-body Hamiltonian within the full configuration-interaction space spanned by 
|to, to') , to = 0, N, to' = N — to. In contrast, the CCSD approach requires here one pa- 
rameter only, c 2 2, which is a solution of a simple algebraic equation of the second degree, see 
Eq. (58). The performance of the CCSD approach for N = 100, 1000 and 10000 interacting 
bosons was tested in terms of three criteria: correlation energy, many-body wavefunction 
and ground-state depletion. It was found that the CCSD is remarkably successful in ob- 
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taining the correlation energy, with absolute error of less than 4% up to a huge coupling 
constant, A = 2 • 10 4 , see Fig. I. The quality of the CCSD many-body wavefunction and its 
ability to accurately describe ground-state depletion were found to be remarkably good for 
all boson numbers and coupling constants as large as A = 10 2 , see Figs. 2-3. For comparison, 
we examined the performance of CISD, which similarly depends on one parameter only, d 2 2- 
CISD was found to be substantially poorer in comparison to CCSD. For instance, it accounts 
for about 50% of the correlation energy only. 

The coupled-cluster theory for bosons presented in this work, as certainly supported by 
the numerical example, is a promising approach to be further developed in the many-boson 
problem. The expressions of the bosonic coupled-cluster theory are much simpler than those 
for fermions since, generally, the ground configuration (permanent) employs one orbital only. 
Consequently, we can treat a very large number of bosons with coupled-cluster expansions 
and employ more virtual (non-occupied) orbitals than the fermionic coupled-cluster can. 
These qualities open the way to study few- to many-boson systems up to a substantial 
interaction where several orbitals are needed to describe the reality. 

The issue of size consistency, as extensively discussed above, is delicate for bosons, and 
depends on the choice of the ground configuration. It relates to the following practical point: 
what is a suitable choice of the ground permanent when a coupled-cluster expansion is to 
be employed with a specific physical system? We can say that, for bosons in a single- well 
trap an useful choice is the simplest permanent where all bosons reside in the same orbital, 
which is the standard mean-field, Gross-Pit aevskii orbital. However, if we wish to usefully 
apply coupled-cluster expansions to a bosonic system undergoing spatial fragmentation or 
superfluid to Mott-insulator transitions, situations that occur in double- and multi-well 
traps, we have to be more careful with the choice of ground configuration, and depart from 
the simplest permanent constructed from the Gross-Pitaevskii mean-field orbital. Recently, 
a more general mean-field theory has been introduced, allowing for bosons to reside in several 
orbitals [27]. We anticipate that in combination with coupled-cluster expansions they can 
be useful for studying many bosons in double- and multi-well traps. 
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APPENDIX A 

The transformed two-body operator V of the transformed Hamiltonian H 
consists of the nine terms listed below. 



V(2 
V(3 
V{A 
V(5 
V(Q 
V(7 
V(8 
V(9 



v = Y. v(p) 

p=l 

M 

E Vnikblblhk 

k=2 
M 

E Vkmblblhh 

k=2 
i M 

- E VnkMhh 

1 k,l=2 
i M 

« E wM&i&i 

Z fc,/=2 
M 

E {Viku + v^blblhk 

k,l=2 
M 

E VijkibWMi 

j,k,l=2 
M 

E tv&K^i 

j,k,l=2 

1 M 

2 E Vmblb^kk 

i,j,k,l=2 
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FIG. 1. (Color online) Performance of CCSD method: correlation energy. Shown is the percent 
of correlation energy, denoted by /^correlation > obtained by the CCSD for N = 100,1000,10000 
bosons and several values of interaction strength Ao- The correlation energy is defined as the differ- 
ence between the Gross-Pitaevskii energy Eqp and the exact energy. The exact energy is obtained 
in our model by diagonalizing the many-body Hamiltonian within the full configuration-interaction 
space. For comparison, the corresponding values obtained by CISD are also plotted. 
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FIG. 2. (Color online) Performance of CCSD method: many-body wavefunction. Shown are 
the coefficients, C2 m , of the normalized many-body wavefunction, J2m=o ^2m \N — 2m, 2m), for 
N = 10000 bosons and A = X (N - 1) = 100 obtained by the CCSD method, see Eqs. (53,54). 
Although the coupling constant is large, it is remarkable that the CCSD C2m coefficients almost 
perfectly match the exact coefficients, namely the red curve "sits" atop the black curve. For 
comparison, the two coefficients of the CISD are also shown, which deviate much from the exact 
solution. 
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FIG. 3. (Color online) Performance of CCSD method: ground-state depletion. Shown is the 
average number of bosons in orbital <p2, {n-zj, for N = 100, 1000, 10000 bosons and several values 
of interaction strength. Because of the different spatial symmetry of <f>\ and cf2, these orbitals are 
the eigenfunctions of the reduced one-particle density matrix (natural orbitals) and the respective 
eigenvalues are (m) and (77-2) with (ni) + (n2) = N. It is seen that the CCSD is extremely successful 
in obtaining the depletion (^2)) which is a very sensitive measure of the exactness of the many-body 
wavefunction, for all ./V up to a large coupling constant, A = 10 2 . The CISD results are also shown 
for comparison. See text for more details. 
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